#include "PowerMethod.h"

#include <iostream>

namespace NumericalAnalysis
{
    // Power method for matrices:
    // The parameters are:
    // Dimension of the matrix (n).
    // File path to read a matrix  (nxn)
    // File path to read a vector (nx1)
    // Max iterations allowed
    // Error tolerance
    // Also you need to pass the eigenvalue and eigenvector 
    // where return the results.
    void PowerMethod(const uint32_t matrixDimension, const std::string& filePathForMatrix, const std::string& filePathForVector, 
    	 const uint32_t maxIterations, const double errorTolerance, double& outEigenvalue,
    	 Utils::Matriz& outEigenvector, uint32_t& outLastIteration)
    {
    	// Build the matrix and vector.
    	Utils::Matriz A = Utils::GetMatrixFromFile(filePathForMatrix, matrixDimension, matrixDimension);
    	Utils::Matriz previousX = Utils::GetMatrixFromFile(filePathForVector, matrixDimension, 1);

    	const double Xp = Utils::GetInfiniteVectorNorm(previousX);
    	previousX = Utils::MultiplyMatrixByLiteral(previousX, 1.0 / Xp);
    	double Yp = 0.0;
    	uint32_t k = 1;
    	Utils::Matriz y(A.mRows, previousX.mColumns);
    	for(; k < maxIterations; ++k)
    	{
			// y = A*x
    		y = Utils::MultiplyMatrices(A, previousX); 
    		
			// Check for zero division.
			// r = ||y|| (Eigenvalue)
			Yp = Utils::GetInfiniteVectorNorm(y);
    		if(Yp == 0.0)
    		{
    			outEigenvalue = Yp;
    			outEigenvector = previousX;
    			
    			break;
    		}

			// x = y / ||y|| (Eigenvector)
    		Utils::Matriz newX = Utils::MultiplyMatrixByLiteral(y, 1.0 / Yp);
    		Utils::Matriz differenceMatrix = Utils::SubtractMatrices(previousX, newX);
    		const double currentError = Utils::GetInfiniteVectorNorm(differenceMatrix);

    		previousX = newX;

    		if(currentError < errorTolerance)
    		{
    			outEigenvalue = Yp;
    			outEigenvector = previousX;
				outLastIteration = k;
    			
    			break;
    		}
    	}

    	outEigenvalue = Yp;
    	outEigenvector = previousX;
    	outLastIteration = k;
    }
}
